Lecture 16

class: center, middle, inverse, title-slide

# Lecture 16
## Factorial Designs
### Psych 10 C
### University of California, Irvine
### 05/06/2022

---

## Factorial designs with 3 factors

- The examples we have worked with have used only two factors or independent
variables, however, some problems will include more than that.

- The steps that we need to take to analyze a `$2\times 2 \times 2$` between 
subjects factorial design are the same as the ones we used for others.

- Let's look at an example about grades on a midterm.

---

## Example

- We want to study the effects of study habits on students' midterm scores.

- We are interested on the effects of three variables:

1. Listening to music while studying.
  
--

2. Distributing the time they spend studying over multiple days vs studying 
  the same amount of time a day before the exam.
  
--

3. Whether students use the provided slides or their own notes.
  
--

- With this 3 factors (music, study and notes) we have a total of 8 independent
groups in the study, for example, we have a group that listened to music, 
studied one day before the exam and used their own notes.

---

## Example

- If we compare all the models that we can build with these 3 factors, we would
end up with a total of 12 different models.

- One Null, 3 main effects models, 4 additive models and 4 interaction models.

- When we have multiple factors in a model comparison setting we have to make 
some decisions. Of course, we can always evaluate every model, but the number 
of models will increase rapidly with the number of factors.

- The second option is to evaluate only the null, main effects and additive 
models.

- The third option (and the better approach) is to have a hypothesis about which
interactions might be relevant to the problem that we are trying to solve and 
only evaluate those models.

---

## Model evaluation

- Given that we only need to compare 8 models, we will just avoid the 
interactions and look at the main effects and additive models.

- The approach will be the same as before. We want to compute the Sum of Squared
Errors, which means that we need the predictions of each model and then the 
error for each observation.

- Given that we are not interested on the full model (interactions) this time we
can start by calculating the grand mean.

```r
null <- scores %>% 
  summarise("pred" = mean(midterm)) %>% 
  pull(pred)
```
--

- The average score in the midterm was equal to 72.4

---

## Main effects models

- Given that we have 3 factors this time (each with two levels), we need to 
estimate and save the effects of each independently, this is because we will
need those values for the additive models.

- First we can obtain the main effects of the **music** factor:

```r
music_mean <- scores %>% 
  group_by(music) %>% 
  summarise("average" = mean(midterm))
alpha <- music_mean$average - null
```

.pull-left[
| music                | average                 |
|:--------------------:|:---------------------:|
| music  | 68.54 |
| no_music | 76.25 |
]
.pull-right[
| Parameter            | Value                 |
|:--------------------:|:---------------------:|
| `$\hat{\alpha}_{music}$`  | -3.85 |
| `$\hat{\alpha}_{no\ music}$` | 3.85 |
]

---

## Main effects models

- Next we will calculate the effects of the **study** factor or how students 
distribute the time they dedicate to study during a week.

```r
study_mean <- scores %>% 
  group_by(study) %>% 
  summarise("average" = mean(midterm))
beta <- study_mean$average - null
```

.pull-left[
| study                | average                 |
|:--------------------:|:---------------------:|
| distributed  | 74.96 |
| one_day | 69.83 |
]
.pull-right[
| Parameter            | Value                 |
|:--------------------:|:---------------------:|
| `$\hat{\beta}_{distributed}$`  | 2.56 |
| `$\hat{\beta}_{one\ day}$` | -2.56 |
]

---

## Main effects models

- The last main effects that we need to calculate are the ones associated to the
type of **notes** that students use to prepare for the midterm.

```r
notes_mean <- scores %>% 
  group_by(notes) %>% 
  summarise("average" = mean(midterm))
gama <- notes_mean$average - null
```

.pull-left[
| notes                | average                 |
|:--------------------:|:---------------------:|
| notes  | 72.89 |
| slides | 71.9 |
]
.pull-right[
| Parameter            | Value                 |
|:--------------------:|:---------------------:|
| `$\hat{\gamma}_{notes}$`  | 0.49 |
| `$\hat{\gamma}_{slides}$` | -0.49 |
]

<br>

- Notice that in comparison to the other effects (with a magnitude between 2 
and 3) the type of notes that students use has a lower effect.

---

## Predictions and errors

- Now that we have the effects of the levels of each of our 3 factors we can 
add the predictions of the model to our data and calculate the errors.

- We will start with the Null and Main effects models.

```r
n_total <- nrow(scores)
scores <- scores %>% 
  mutate("prediction_null" = null,
         "error_null" = (midterm - null)^2)
sse_null <- sum(scores$error_null)
mse_null <- 1/n_total * sse_null
bic_null <- n_total * log(mse_null) + 1 * log(n_total)
```

- The SSE of the null model was equal to 5658.44.

---

## Predictions and error main effects

- The predictions of the **music** main effects model are: 
`$$\hat{\mu} + \hat{\alpha}_{music} \quad \hat{\mu} + \hat{\alpha}_{no\ music}$$`

```r
scores <- scores %>% 
  mutate("prediction_music" = case_when(music == "music" ~ null + alpha[1],
                                        music == "no_music" ~ null + alpha[2]),
         "error_music" = (midterm - prediction_music)^2)
sse_music <- sum(scores$error_music)
mse_music <- 1/n_total * sse_music
r2_music <- (sse_null - sse_music)/sse_null
bic_music <- n_total * log(mse_music) + 2 * log(n_total)
```

- The SSE for the **music** main effects model was 3519.38

---

## Predictions and error main effects

- Now we will calculate the predictions of the **study** main effects model, 
this is the effect of distributing study time across multiple days before a 
test or if they study one day before the test.

- The predictions of the **study** main effects model are:
`$$\hat{\mu} + \hat{\beta}_{distributed} \quad \hat{\mu} + \hat{\beta}_{one\ day}$$`
--

```r
scores <- scores %>% 
  mutate("prediction_study" = case_when(study == "distributed" ~ null + beta[1],
                                        study == "one_day" ~ null + beta[2]),
         "error_study" = (midterm - prediction_study)^2)
sse_study <- sum(scores$error_study)
mse_study <- 1/n_total * sse_study
r2_study <- (sse_null - sse_study)/sse_null
bic_study <- n_total * log(mse_study) + 2 * log(n_total)
```

- The SSE of the **study** main effects model was 4712.88

---

## Predictions and error main effects

- The last main effects model is the main effects of **notes** or whether 
students used their own notes or the notes provided by the instructor.

- The predictions of the **notes** model are:
`$$\hat{\mu} + \hat{\gamma}_{slides} \quad \hat{\mu} + \hat{\gamma}_{notes}$$`

```r
scores <- scores %>% 
  mutate("prediction_notes" = case_when(notes == "notes" ~ null + gama[1],
                                        notes == "slides" ~ null + gama[2]),
         "error_notes" = (midterm - prediction_notes)^2)
sse_notes <- sum(scores$error_notes)
mse_notes <- 1/n_total * sse_notes
r2_notes <- (sse_null - sse_notes)/sse_null
bic_notes <- n_total * log(mse_notes) + 2 * log(n_total)
```

- The SSE for the **notes** main effects model was 5623.43

---

## Model evaluation

- Before we look at the results from the additive models, we can look at the
comparison between the main effects models and the null.

- Again, we can summarize our results using a table.

- From the table we can see that the best model for now is the main effects of 
**music** given that it has the lowest BIC.

- However, an important thing to note is that the BIC of the **notes** main 
effects model is larger than the Null.

---

## Model evaluation

- The fact that the **notes** main effects model performs worse than the Null
indicates that, knowing whether a student used the slides provided or 
their own notes to study, does not improve the accuracy of our predictions
in comparison to using the mean of all participants as a prediction (the grand 
mean).

- Given that additive models assume that the effects of each factor are 
independent it is very likely that models that include the **notes** factor will
not perform better than other additive models.

- For example, it is likely that the **music** main effects model is going to 
be better (will have a lower BIC) than the additive model that assumes 
that both **music** and **notes** have an effect on midterm scores (the additive 
model that only includes those two factors).

---

## Predictions and errors, additive models

- As we said before, we have 4 additive models that we need to evaluate in a 
`$2\times 2 \times 2$` between subjects factorial design.

- The first model assumes that only factor 1 and factor 2 have an effect on the 
expected value of the dependent variable.

- In our example, this would be the additive model that includes the factor 
**music** and **study**.

- We already have the main effects of each ( `$\alpha$` and `$\beta$` respectively)
so we can add the predictions of the model to the data.

```r
scores <- scores %>% 
  mutate("prediction_additive_ms" = case_when(
    music == "no_music" & study == "distributed" ~ null + alpha[2] + beta[1],
    music == "no_music" & study == "one_day" ~ null + alpha[2] + beta[2],
    music == "music" & study == "distributed" ~ null + alpha[1] + beta[1],
    music == "music" & study == "one_day" ~ null + alpha[1] + beta[2]),
         "error_additive_ms" = (midterm - prediction_additive_ms)^2)
```

---

## Music and study additive model

- Now with the predictions of the model we can calculate the SSE and the other
values that we need for the comparison:

```r
sse_ms <- sum(scores$error_additive_ms)
mse_ms <- 1/n_total * sse_ms
r2_ms <- (sse_null - sse_ms)/sse_null
bic_ms <- n_total * log(mse_ms) + 3 * log(n_total)
```

- The SSE of the additive model that includes the **music** and **study** 
factors was 2573.81

---

## Music and notes additive model

- The second additive model assumes that **only** the first and third factor 
have an effect on the expected value of our dependent variable.

- In our example, the model assumes that only **music** and **notes** have an 
effect on the expected midterm score of a student. We can add the predictions
of the model to our data using a similar code as before:

```r
scores <- scores %>% 
  mutate("prediction_additive_mn" = case_when(
    music == "no_music" & notes == "notes" ~ null + alpha[2] + gama[1],
    music == "no_music" & notes == "slides" ~ null + alpha[2] + gama[2],
    music == "music" & notes == "notes" ~ null + alpha[1] + gama[1],
    music == "music" & notes == "slides" ~ null + alpha[1] + gama[2]),
         "error_additive_mn" = (midterm - prediction_additive_mn)^2)
sse_mn <- sum(scores$error_additive_mn)
mse_mn <- 1/n_total * sse_mn
r2_mn <- (sse_null - sse_mn)/sse_null
bic_mn <- n_total * log(mse_mn) + 3 * log(n_total)
```

- The SSE of the additive model that includes the **music** and **notes** 
factors was 3484.37

---

## Study and notes additive model

- The last additive model that includes **only** two of the factors is the model
that assumes that only factor 2 and factor 3 have an effect on the expected 
value of the dependent variable.

- In our example, the model assumes that only the factors **study** and **notes**
have an effect on the expected midterm score of a student.

```r
scores <- scores %>% 
  mutate("prediction_additive_sn" = case_when(
    study == "distributed" & notes == "notes" ~ null + beta[1] + gama[1],
    study == "distributed" & notes == "slides" ~ null + beta[1] + gama[2],
    study == "one_day" & notes == "notes" ~ null + beta[2] + gama[1],
    study == "one_day" & notes == "slides" ~ null + beta[2] + gama[2]),
         "error_additive_sn" = (midterm - prediction_additive_sn)^2)
sse_sn <- sum(scores$error_additive_sn)
mse_sn <- 1/n_total * sse_sn
r2_sn <- (sse_null - sse_sn)/sse_null
bic_sn <- n_total * log(mse_sn) + 3 * log(n_total)
```

- The SSE of the additive model that includes the **study** and **notes** 
factors was 4677.87

---

## Additive model with 3 factors

- The last model that we need to consider (if we ignore the interactions) is the 
additive model that assumes that all 3 factors have an effect on the expected
value of our dependent variable.

- The predictions of this model now will include the estimate of the grand mean 
`$\hat{\mu}$`, the main effect of **music** `$\hat{\alpha}$`, the main effect of 
**study** `$\hat{\beta}$` and the main effect of **notes** `$\hat{\gamma}$`.

- Now we need to check 3 conditions instead of 2, this is because a participant 
could have been exposed to a specific combination of the levels of the 3 factors. 
For example, for the combination **no_music**-**distributed**-**slides** the 
prediction should be:
`$$\hat{\mu} + \hat{\alpha}_{no\ music} + \hat{\beta}_{distributed} + \hat{\gamma}_{slides}$$`

- However, the prediction for a student that was exposed to the levels **music**-
**distributed**-**slides** would be: 
`$$\hat{\mu} + \hat{\alpha}_{music} + \hat{\beta}_{distributed} + \hat{\gamma}_{slides}$$`

---

### Additive model with 3 factors predictions and error

```r
scores <- scores %>% 
  mutate("prediction_additive" = case_when(
    music == "music" & study == "distributed" & notes == "notes" ~ 
      null + alpha[1] + beta[1] + gama[1],
    music == "music" & study == "distributed" & notes == "slides" ~ 
      null + alpha[1] + beta[1] + gama[2],
    music == "music" & study == "one_day" & notes == "notes" ~ 
      null + alpha[1] + beta[2] + gama[1],
    music == "music" & study == "one_day" & notes == "slides" ~ 
      null + alpha[1] + beta[2] + gama[2],
    music == "no_music" & study == "distributed" & notes == "notes" ~ 
      null + alpha[2] + beta[1] + gama[1],
    music == "no_music" & study == "distributed" & notes == "slides" ~ 
      null + alpha[2] + beta[1] + gama[2],
    music == "no_music" & study == "one_day" & notes == "notes" ~ 
      null + alpha[2] + beta[2] + gama[1],
    music == "no_music" & study == "one_day" & notes == "slides" ~ 
      null + alpha[2] + beta[2] + gama[2]),
         "error_additive" = (midterm - prediction_additive)^2)
sse_additive <- sum(scores$error_additive)
mse_additive <- 1/n_total * sse_additive
r2_additive <- (sse_null - sse_additive)/sse_null
bic_additive <- n_total * log(mse_additive) + 4 * log(n_total)
```

---

## Comparing additive models

- Now that we have all the values we need we can add the additive models to 
the table that we started before:

| Model | Parameters | MSE | `$R^2$` | BIC |
|-------|:----------:|:---:|:-----:|:---:|
| Null | 1 | 39.29 | NA | 533.61 |
| ME Music | 2 | 24.44 | 0.38 | 470.2 |
| ME Study | 2 | 32.73 | 0.17 | 512.25 |
| ME Notes | 2 | 39.05 | 0.01 | 537.68 |
| Music + Study | 3 | 17.87 | 0.55 | 430.11 |
| Music + Notes | 3 | 24.2 | 0.38 | 473.73 |
| Study + Notes | 3 | 32.49 | 0.17 | 516.14 |
| Music + Study + Notes | 4 | 17.63 | 0.55 | 433.11 |

- Like we expected, adding the **notes** factor does not improve on the models
that where already better at accounting for the data, for example, the model that 
includes **study** and **notes** has a higher BIC than the model that only 
includes **study**